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Abstract 



This paper contains a numerical stability analysis of factorization algorithms for com- 
puting the Cholesky decomposition of symmetric positive definite matrices of displacement 
rank 2. The algorithms in the class can be expressed as sequences of elementary downdating 
steps. The stability of the factorization algorithms follows directly from the numerical prop- 
erties of algorithms for realizing elementary downdating operations. It is shown that the 
Bareiss algorithm for factorizing a symmetric positive definite Toeplitz matrix is in the class 
and hence the Bareiss algorithm is stable. Some numerical experiments that compare behav- 
ior of the Bareiss algorithm and the Levinson algorithm are presented. These experiments 
indicate that in general (when the reflection coefficients are not all positive) the Levinson 
algorithm can give much larger residuals than the Bareiss algorithm. 



1 Introduction 

We consider the numerical stability of algorithms for solving a linear system 

Tx = b, (1.1) 

where T is an n x n positive definite Toeplitz matrix and b is an n x 1 vector. We assume that 
the system is solved in floating point arithmetic with relative precision e by first computing the 
Cholesky factor of T. Hence the emphasis of the paper is on factorization algorithms for the 
matrix T. 

Roundoff error analyses of Toeplitz systems solvers have been given by Cybenko [10] and 
Sweet [22]. Cybenko showed that the Levinson-Durbin algorithm produces a residual which, 
under the condition that all reflection coefficients are positive, is of comparable size to that 
produced by the well behaved Cholesky method. He hypothesised that the same is true even 
if the reflection coefficients are not all positive. If correct, this would indicate that numerical 
quality of the Levinson-Durbin algorithm is comparable to that of the Cholesky method. 
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In his PhD thesis [22], Sweet presented a roundoff error analysis of a variant of the Bareiss 
algorithm [2], and concluded that the algorithm is numerically stable (in the sense specified in 
Section 7). In this paper we strengthen and generalize these early results on the stability of 
the Bareiss algorithm. In particular, our approach via elementary downdating greatly simplifies 
roundoff error analysis and makes it applicable to a larger-than-Toeplitz class of matrices. 

After introducing the notation and the concept of elementary downdating in Sections 2 and 3, 
in Section 4 we derive matrix factorization algorithms as a sequence of elementary downdating 
operations (see also [4]). In Section 5 we present a first order analysis by bounding the first 
term in an asymptotic expansion for the error in powers of e. By analyzing the propagation of 
first order error in the sequence of downdatings that define the algorithms, we obtain bounds 
on the perturbations of the factors in the decompositions. We show that the computed upper 
triangular factor U of a positive definite Toeplitz matrix T satisfies 

T = U T U + AT , ||AT|| < c{n)e\\T\\ , 

where c(n) is a low order polynomial in n and is independent of the condition number of T. 
Many of the results of Sections 2-5 were first reported in [5] , which also contains some results 
on the stability of Levinson's algorithm. 

In Section 6 we discuss the connection with the Bareiss algorithm and conclude that the 
Bareiss algorithm is stable for the class of symmetric positive definite matrices. Finally, in 
Section 7 we report some interesting numerical examples that contrast the behaviour of the 
Bareiss algorithm with that of the Levinson algorithm. We show numerically that, in cases 
where the reflection coefficients are not all positive, the Levinson algorithm can give much 
larger residuals than the Bareiss or Cholesky algorithms. 



2 Notation 



Unless it is clear from the context, all vectors are real and of dimension n. Likewise, all matrices 
are real and their default dimension is n x n. If a £ !R n , ||a|| denotes the usual Euclidean norm, 
and if T £ K nxri , ||T|| denotes the induced matrix norm: 

IITII = max llTall . 

Ha||=l" 
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Our primary interest is in a symmetric positive definite Toeplitz matrix T whose i, jth entry 

Uj = t\i-j\- 

We denote by e^, k = 1, . . . ,n, the unit vector whose kth. element is 1 and whose other 
elements are 0. We use the following special matrices: 



n-1 

z = e *+ief 
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dn—i+l^i 
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The matrix Z is known as a shift-down matrix. We also make use of powers of the matrix Z, 
for which we introduce the following notation: 



I / iffc = 0, 
k ~ \ Z k if k > 0. 

The antidiagonal matrix J is called a reversal matrix, because the effect of applying J to a 
vector is to reverse the order of components of the vector: 



J 









X2 




•En— 1 









The hyperbolic rotation matrix H(9) £ K 2x2 is defined by 

1 



H{9) 

The matrix H (9) satisfies the relation 

H(9) 



cos 9 



1 — sin ( 
sin0 1 



(2.1) 



1 
-1 



H(9) 



1 
-1 



and it has eigenvalues \\{9), \2{9) given by 

\ 1 (9) = \ 2 1 (9) = secfl-tan#. 



(2.2) 



For a given pair of real numbers a and b with |a| > \b\, there exists a hyperbolic rotation matrix 
H (9) such that 



H{0) 



The angle of rotation 9 is determined by 



a 




vV - 6 2 


b 








. b , Va^W 

sin 6> = - , cos = 

a a 



(2.3) 



(2.4) 



3 Elementary Downdating 

In this section we introduce the concept of elementary downdating. The elementary downdating 
problem is a special case of a more general downdating problem that arises in Cholesky factor- 
ization of a positive definite difference of two outer product matrices [1, 6, 7, 12]. In Section 4, 
factorization algorithms are derived in terms of a sequence of downdating steps. The numerical 
properties of the algorithms are then related to the properties of the sequence of elementary 
downdating steps. 

Let Ufc, Vfc G ffl 1 have the following form: 

k 
i 

ul = [0 ... x 
v£ = [0 ... 



X X ... X ] , 
X X ... X ] , 

t 

k + 1 
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that is: 

ejufc = , j < k , and 
Applying the shift-down matrix Z to u^, we have 



T 



[0 
[0 
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k + 1 

I 
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X 

t 

k + 1 



0,j<k 
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X 



xj 
xl 



Suppose that we wish to find u^+i, v k+ i G W 1 to satisfy 



T T T T 

u fc+i u fc+i - v fc+1 v fc+1 = Zu k u k Z 



(3.1) 



where 



T 

u fc+i 

T 



[0 
[0 








k + 1 

I 
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x 
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t 

+ 2 



xj 
xl 



that is 



e ■ u fe+ i = , j < k + 1 , and e ■ v fc+i = , j < k + 1 



(3.2) 



We refer to the problem of finding u k+ i and v k+ i to satisfy (3.1), given and v^, as the 
elementary downdating problem. It can be rewritten as follows: 



[ufc+i v fc+i] 



1 
-1 



T 
U k+1 

T 
v fc+l 



[Zu k v fc ] 



1 

-1 



u; 



From (2.1), (2.3) and (2.4), it is clear that the vectors u k+ i and v k+i can be found by using a 
hyperbolic rotation H (9 k ) defined by the following relations: 



sin6> fc = e fe+1 v fe /e fc u fc , 
cos 6 k = \jl — sin 2 9k , 



and 



T 
u fe+l 

T 
V fc+1 



(3.3a) 
(3.3b) 

(3.4) 



■H{6 k ) 

The elementary downdating problem has a unique solution (up to obvious sign changes) if 

\elu k \ > \el +1 v k \ . 

The calculation of u k+ \, v k+ \ via (3.4) can be performed in the obvious manner. Following 
common usage, algorithms which perform downdating in this manner will be referred to as 
hyperbolic downdating algorithms. 

Some computational advantages may be obtained by rewriting (3.1) as follows: 



u 



fc+1 



[Zu k v fe+ i] 



T 
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Consider now an orthogonal rotation matrix G(9k), 



G{9 k 



cos 6k sin 9 k 
- sin Ok cos 9k 



where cos9k and sin^ are defined by (3.3b) and (3.3a), respectively. Then it is easy to check 
that 



or, equivalently, 



Thus, we may rewrite (3.6) as 



G(9 k ) 


T 




T 










T 
u k+l 

T 


= G(9k) T 


" <z T ' 

T 











Vfe+i = (v fe - sm9 k Zu k )/ cos 9 k , 
u fc+ i = -sin6> fc v fe+ i +cosOkZu k • 



(3.5) 



(3.6) 



(3.7a) 
(3.7b) 



Note that equation (3.7a) is the same as the second component of (3.4). However, (3.7b) differs 
from the first component of (3.4) as it uses Vk+i in place of v& to define Uk+i- It is possible to 
construct an alternative algorithm by using the first component of (3.5) to define U&+1. This 
leads to the following formulas: 



Ufe+i = (Zu k - sin 9 k v k )/ cos 9 k , 
Vfc+i = -sin9 k u k+l +cos(9 fc v fc . 



(3.8a) 
(3.8b) 



We call algorithms based on (3.7a)-(3.7b) or (3.8a)-(3.8b) mixed elementary downdating al- 
gorithms. The reason for considering mixed algorithms is that they have superior stability 
properties to hyperbolic algorithms in the following sense. 

Let iifc, Vfc be the values of u^, that are computed in floating point arithmetic with relative 
machine precision e. The computed values u^, satisfy a perturbed version of (3.1), that is, 



Ufc+iu fe+1 - v fc+1 v fc+1 = Zu k u k Z 



v fe v^ + eG k + Q(e 2 



(3.9) 



where the second order term 0(e 2 ) should be understood as a matrix whose elements are bounded 
by a constant multiple of e 2 ||Gfc||. The norm of the perturbation Gk depends on the precise 
specification of the algorithm used. It can be shown [6] that the term Gk satisfies 



\Gk\\ < c ri 



Zu k \\ 2 + ||v fc || 2 + ||u fc+ i|| 2 + ||v fc+ i|| 2 ) 



(3.10) 



when a mixed downdating strategy is used (here c m is a positive constant). When hyperbolic 
downdating is used the term Gk satisfies 



||Gfc|| < c h ||ir(0 fc )|| (||Zufc|| + ||vfc||) (||u fc+ i|| + ||v fc+ i||) , (3.11) 

where Ch is a positive constant [6]. (The constants c m and Ch are dependent on implementation 
details, but are of order unity and independent of n.) Note the presence of the multiplier \\H(9k) \\ 
in the bound (3.11) but not in (3.10). In view of (2.2), ||i?(^)|| could be large. The significance 
of the multiplier ||i?(^)|| depends on the context in which the downdating arises. We consider 
the implications of the bounds (3.10) and (3.11) in Section 5 after we make a connection between 
downdating and the factorization of Toeplitz matrices. 
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It is easily seen that a single step of the hyperbolic or mixed downdating algorithm requires 
4(n— k)+0(l) multiplications. A substantial increase in efficiency can be achieved by considering 
the following modified downdating problem. Given aj, 6 3f and w&, x k G 3? n that satisfy 

ejw fc = , j < k and ejx fc = , j < k , 

find afc+i, Pk+i and w^+i, x& + i £ !R n that satisfy 

«I+i w fc+i w fc+i - /5i+i x fe+i x fc+i = a 2 k Zw k wlZ T - f3lx k xl , 

with 

e J w fc = , j < k and ejx^. = , j < k . 
If we make the identification 



u fc = Q fc w fc and v k = f3 k x k , 



then we find that the modified elementary downdating problem is equivalent to the elementary 
downdating problem. However, the extra parameters can be chosen judiciously to eliminate 
some multiplications. For example, if we take a k = f3 k , a k+ \ = fik+i, then from (3.3a), (3.3b) 
and (3.4), 



and 



sin6» fe = el +1 x k /elw k , 
afc+i = <W cos9 k , 

w fc+ i = Zw k - sin 6 k x k , 
Xfc+i = - sin 9 k Zw k + x fc . 



(3.12a) 
(3.12b) 



(3.13a) 
(3.13b) 



Equations (3.12a)-(3.13b) form a basis for a scaled hyperbolic elementary downdating algorithm 
which requires 2(n — A;) + 0(1) multiplications. This is about half the number required by the 
unsealed algorithm based on (3.4). (The price is an increased likelihood of underflow or overflow, 
but this can be avoided if suitable precautions are taken in the code.) 

Similarly, from (3.7a) and (3.7b) we can obtain a scaled mixed elementary downdating algo- 
rithm via 

sin6» fc = /3 k el +1 x k /a k elw k , 
"fe+i = a k cos6 k , 
(3 k+1 = f3 k / cos 9 k , 



and 

sin6 k a k 
x fc+i = x fc p Zw k , 

Pk 

sm9 k f3 k+ i 
Wfe+i = x fe+ i + Zw fc . 

The stability properties of scaled mixed algorithms are similar to those of the corresponding 
unsealed algorithms [12]. 
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4 Symmetric Factorization 



We adopt the following definition from [18]. 

Definition 4.1: An n x n symmetric matrix T has displacement rank 2 iff there exist vectors 
u, v G W 1 such that 

T - ZTZ T = uu T - vv T . (4.1) 

□ 

The vectors u and v are called the generators of T and determine the matrix T uniquely. 
Whenever we want to stress the dependence of T on u and v we write T = T(u , v). 

In the sequel we will be concerned with a subset T of all matrices satisfying (4.1). The 
subset is defined as follows. 
Definition 4.2: A matrix T is in T iff 

(a) T is positive definite, 

(b) T satisfies (4.1) with generators u and v, 

(c) v T ei = 0, i.e., the first component of v is zero. 



□ 



It is well known that positive definite n x n Toeplitz matrices form a subset of T ■ Indeed, if 
T = (t li _ jl )U^, then 



where 



u 



v 



T - ZTZ T = uu T - vv T 



— (to , h , ■ ■ ■ , tn-l) I Vto , 

= (0 , t\ , . . . , t„_i) /Vto . 



The set T also contains matrices which are not Toeplitz, as the following example shows. 
Example: Let 





25 


20 


15 




5 







T = 


20 


32 


29 


. u = 


4 


and v = 


3 




15 


29 


40 




3 




1 



It is easy to check that T is positive definite. Moreover, 





' 25 


20 


15 " 




' 25 


20 


15 " 




" 
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T - ZTZ T = 


20 
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9 




20 


16 


12 
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3 




15 
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8 




15 


12 


9 







3 


1 



T T 
uu — vv 



Hence T = T(u , v) <E T, but T is not Toeplitz. 

□ 

We now establish a connection between the elementary downdating problem and symmetric 
factorizations of a matrix from the set T ■ 
Let T = T(u , v) G T. Set 

Ui = U, Vi = V 

and, for k = 1, . . . , n — 1, solve the elementary downdating problem defined by (3.1), 

u fc+ iu£ +1 - v fe+1 v^ +1 = Zu k ulZ T - w k wl , 
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which we assume for the moment has a solution for each k. On summing over k = 1, . . . , n — 1 
we obtain 

n— 1 n— 1 n— 1 ra— 1 

X! u fc+i u fe+i - H v fc+i v fe+i = £ Zu k ulZ T - v fe v^ . 
fc=l fc=l fc=l fc=l 

If we now observe that, from (3.2), 

Zu n = v n = , 

we arrive at the following relation: 

n / n \ 

£ u fc u£ - Z £ u fe u£ Z T = Ul uf - v lV f , (4.2) 

fc=i \fc=i / 

which implies that J2k=i u k u k e Moreover, as matrices having the same generators are 
identical, we obtain 

n 

T = ]T u fe u£ = U T U , 

where 

n 

U = J2 e k^k 

k=l 

is upper triangular, and hence is the Cholesky factor of T. We have derived, albeit in a rather 
indirect manner, the basis of an algorithm for calculating the Cholesky decomposition of a matrix 
from the set T. 

We now return to the question of existence of a solution to the elementary downdating 
problem for each k = 1, . . . ,n — 1. It is easy to verify that, if T G T, then \ej ui| > le^vij. 
Using (4.2) and (3.1), it can be shown by induction on k that 

\elu k \ > \el +1 v k \, k = 2,...,n — l. 

Consequently, | sin#fc| < 1 in (3.3a), and the elementary downdating problem has a solution for 
each k = 1, . . . , n — 1. 

To summarize, we have the following algorithm for factorizing a matrix T = T(u , v) G T ■ 
Algorithm FACTOR(T): 
Set ui = u, vi = v. 

For k = 1, . . . , n — 1 calculate u k+ i, v^ +1 such that 



u fc+i u fc+i - V fc+l V fc+l = ^UfcU fc Z -v fe v fc , 
e fe+i v fc+i = . 

Then T = U T U, where U = £Li e fc u£. 

□ 

In fact we have not one algorithm but a class of factorization algorithms, where each al- 
gorithm corresponds to a particular way of realizing the elementary downdating steps. For 
example, the connection with the scaled elementary downdating problem is straightforward. On 
making the identification 

u k = a k w k and v fe = /3 fc x fc , (4.3) 



we obtain 



T = W T D 2 W , 
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where 



n 

w = Y, e ^k, 

fc=i 

n 

D = Y «fc e fc e fc • 
fc=l 

It is clear from Section 3 that Algorithm FACTOR(T) requires 2n 2 + 0(n) multiplications 
when the unsealed version of elementary downdating is used, and n 2 + 0{n) multiplications 
when the scaled version of elementary downdating is used. However, in the sequel we do not 
dwell on the precise details of algorithms. Using (4.3), we can relate algorithms based on the 
scaled elementary downdating problem to those based on the unsealed elementary downdating 
problem. Thus, for simplicity, we consider only the unsealed elementary downdating algorithms. 

5 Analysis of Factorization Algorithms 

In this section we present a numerical stability analysis of the factorization of T £ T via 
Algorithm FACTOR(T). The result of the analysis is applied to the case when the matrix T is 
Toeplitz. 

Let Ufc, Vfc be the values of u&, that are computed in floating point arithmetic with relative 
machine relative precision e. The computed quantities and v^ satisfy the relations 

ti fc = u fc + 0(e), Vfc = v fc + 0(e), (5.1) 

and the aim of this section is to provide a first order analysis of the error. By a first order analysis 
we mean that the error can be bounded by a function which has an asymptotic expansion in 
powers of e, but we only consider the first term of this asymptotic expansion. One should think 
of e — > 0+ while the problem remains fixed [19]. Thus, in this section (except for Corollary 5.1) 
we omit functions of n from the "O" terms in relations such as (5.1) and (5.2). 

The computed vectors Ufc, v& satisfy a perturbed version (3.9) of (3.1). On summing (3.9) 
over k = 1, . . . , n — 1 we obtain 

n— 1 

f - ZTZ T = Qiuf - vivf - (2u„u^ T - v n v£) + e ]T G k + 0(e 2 ) , 

k=i 

where 

f = U T U , 

n 

U = e ^k ■ 

k=l 

Since 

Zu n = 0(e), v n = O(e) , 

we find that 

n-l 

f - ZTZ T = uiuf - vi vf + e G k + 0(e 2 ) . (5.2) 

fc=i 

Now define 

E = f — T. (5.3) 
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Then, using (4.1), (5.2) and (5.3), 



n— 1 



E - ZEZ T = uiuf - uu T + vivf - vv T + e ^ G fe + 0(e 2 ) . 



fc=i 



In a similar manner we obtain expressions for ZjEZj — Zj +1 EZj +l , j = 0, . . . , n — 1. Summing 
over j gives 

n— 1 n— 1 n— 1 

E=Y,Z j ((uiuf - uiuf ) + (vivf - v lV f )) + e E E ^' G ^J + 0(e 2 ) . (5.4) 

j=0 j=0 fe=l 

We see from (5.4) that the error consists of two parts - the first associated with initial errors 
and the second associated with the fact that (5.2) contains an inhomogeneous term. Now 



||uiuf - uu T || < 2||u|| ||tii - u|| + 0(e 2 ) , 
||vivf - vv T || < 2 1 1 v 1 1 || vi - v|| + 0(e 2 ) . 

Furthermore, from (4.1), 

Tr(T) - Tr(ZTZ T ) = ||u|| 2 - ||v|| 2 > , 

and hence 
n-l 

|E^(uiur-uu T + Vl vf - vv T )Zj|| <2n||u||(||ui-u|| + ||vi-v||) + 0(e 2 ) . (5.5) 

3=0 

This demonstrates that initial errors do not propagate unduly. To investigate the double sum 
in (5.4) we require a preliminary result. 

Lemma 5.1 For k = 1, 2, ... , n — 1 and j = 0, 1, 2, . . . , 

||ZjVfc|| < H-Zj+iUfcH . 



□ 



Proof Let 



It is easy to verify that 



k n 

T k = T-j2"i"T = E 

1=1 l=k+l 



T k - ZT k Z T = Zu k ulZ T - w k wl 



and, since T k is positive semi-definite, 

Tr{zp k Z] - Z 3+1 T k Zj +1 ) = H^+iUfcll 2 - ||Z iVfc || 2 > . 

□ 

We now demonstrate stability when the mixed version of elementary downdating is used in 
Algorithm FACTOR(T). In this case the inhomogeneous term G k satisfies a shifted version 
of (3.10), that is 

\\ZjG k Zj\\ < c m (||Z i+1 u fc || 2 + ll^-vfcll 2 + ll^-ufe+ill 2 + ll^-vfc+ill 2 ) , (5.6) 
where c m is a positive constant. 
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Theorem 5.1 Assume that (3.9) and (5.6) hold. Then 



n-l 



\T - U T U\\ < 2n\\u\\ (||£ii - u|| + ||vi - v||) + 4ec m ^ Tr(ZjTZj) + 0(e 2 ) 



3=0 

Proof Using Lemma 5.1, 



□ 



Furthermore, since 



it follows that 



\ZjGkZj\\ < 2c m (\\Zj + iu.k\\ 2 + \\ZjUk+i\\ 2 ^J . 
Tr{Z 3 TZj) = J2 \\Z 3 u k \\ 2 , 



k=i 

n—l n n—l 



\zZzZ Z 3GkZj <Ac m Y,Tr(Z 3 TZj). (5.7) 

3=0 k=l j=0 

The result now follows from (5.4), (5.5) and (5.7). 

□ 

For the hyperbolic version of the elementary downdating algorithms a shifted version of the 
weaker bound (3.11) on G k holds (see [6]), namely 

\\ZjGkZfW < c h \\H(e k )\\(\\Z j+1 u k \\ + IIZjVfcllXHZjUfc+iH + ||Z,-v fc+ i||) . (5.8) 

By Lemma 5.1, this simplifies to 

\\Z 3 G k Zj\\ < 4^11^(^)11 HZj+iUfcH \\ZjU k+1 \\ . (5.9) 

The essential difference between (3.10) and (3.11) is the occurence of the multiplier ||iJ(^)|| 
which can be quite large. This term explains numerical difficulties in applications such as the 
downdating of a Cholesky decomposition [6]. However, because of the special structure of the 
matrix T, it is of lesser importance here, in view of the following result. 
Lemma 5.2 For k = 1, 2, ... , n — l and j = 0, 1, . . . , n — k, 



\H(9 k )\\ \\Z jUk+1 \\ <2(n-k-j)\\Z J+1 u k \\. 



□ 



Proof It is easy to verify from (3.4) that 

1 T sin 6 k 



cos 9 k 

and from (2.1) that 

\\H{9 k )\\ = 



(Ujfc+i T v fe+ i) = Zu k T v fe 
1 + |sin0| 



cos 6 
Thus, 

\\H(9 k )\\ ||ZjU fc+ i|| < \\H(9 k )\\ HZjVfc+iH + ll^+iUfeU + HZjVfcH 
< \\H(9 k )\\ ||Z,- +1 u fc+ i|| + 2||Z i ,- +1 u fc || , 

where the last inequality was obtained using Lemma 5.1. Thus 

n— k 

\\H{9 k )\\ ||Z,-u fc+ i|| < 2 J2 \\Zi"k\\ , 

l=j+l 

and the result follows. □ 
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Remark Lemma 5.2 does not hold for the computed quantities unless we introduce an O(e) 
term. However in a first order analysis we only need it to hold for the exact quantities. 

Theorem 5.2 Assume that (3.9) and (5.8) hold. Then 

n-l 

\\T - U T U\\ < 2n||u||(||ui - u|| + || Vl - v||) + 8ec^(n - j)Tr{Z 3 TZj) + 0(e 2 ) . 

3=1 



□ 



Proof Applying Lemma 5.2 to (5.9) gives 

\\ZjG k Zj\\ < 8c h (n - j - l)||^+iu fc || 2 , 



and hence 



n-ln-l n-ln-1 

\j2zZ Z 3 G kZf\\ < 8c h J2J2( n -3)\\Z J u k \\ 2 

j=0 k=l 3=1 k=l 

n-1 

< 8c h J2(n-j)Tr(ZjTZj). (5.10) 

3=1 



□ 



The result now follows from (5.4), (5.5) and (5.10). 
Note that, when T is Toeplitz, 

Tr(Z 3 TZj) = {n- j)t . 

Hence, from Theorems 5.1 and 5.2, we obtain our main result on the stability of the factorization 
algorithms based on Algorithm FACTOR(T) for a symmetric positive definite Toeplitz matrix: 
Corollary 5.1 The factorization algorithm FACTOR(T) applied to a symmetric positive 
definite Toeplitz matrix T produces an upper triangular matrix U such that 

T = U T U + AT , 

where || AT|| = O(eton 2 ) when mixed downdating is used, and || AT|| = O(eton 3 ) when hyperbolic 
downdating is used. 

□ 



6 The Connection with the Bareiss algorithm 

In his 1969 paper [2], Bareiss proposed an 0(n 2 ) algorithm for solving Toeplitz linear systems. 
For a symmetric Toeplitz matrix T, the algorithm, called a symmetric Bareiss algorithm in [22], 
can be expressed as follows. Start with a matrix := T and partition it in two ways: 

where is the first row of T and is the last row of T. Now, starting from A^°\ compute 
successively two matrix sequences {^W} and {A(~^}, i = 1, . . . ,n— 1, according to the relations 

A (i) = A (i-i) _ ai _ lZi A(- i+1 ^ , = A^ + ^ - a„ i+1 ZfA^ . (6.1) 
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For 1 < i < n — 1, partition A.W and l ) as follows: 

AW - f ^ ^ A™ - ( T( ~ i_1) 



where £/W denotes the first i + 1 rows of A^ , and denotes the last i + 1 rows of A*"*). It 
is shown in [2] that 

(a) T( i+1 ) and T^" 1 ) are Toeplitz, 

(b) for a proper choice of CKj-i and a-j+i, the matrices LW and f/W a re lower and upper 

trapezoidal, respectively, and 

(c) with the choice of aj_i and a_j+i as in (b), the Toeplitz matrix T^ 1 ^ 1 ^ has zero elements 

in positions 2, . . . , i + 1 of its first row, while the Toeplitz matrix T^ l+1 ^ has zero elements 
in positions n — 1, . . . , n — i of its last row. 



Pictorially, 



/ x 




T (i+i) 



X 



X X 

x 



A(- 



T (-i-i) 



x 

v x ••• 
/ x ( 



X 
X 







x \ 



V x 




x / 



After n — 1 steps, the matrices A( n_1 ) and A(~ ra+1 ) are lower and upper triangular, respec- 
tively. At step i only rows i + 1, . . . , n of A^ and rows 1, 2, . . . , n — i of A^~ 1 ^ are modified; the 
remaining rows stay unchanged. Moreover, Bareiss [2] noticed that, because of the symmetry 
of T, 

= J i+1 T ( - l -Vj n and = a_ i+1 , (6.2) 

Here Jj+i and J n are the reversal matrices of dimension (i + 1) x {i + 1) and n x n respectively. 

Now, taking into account (6.2), it can be seen that the essential part of a step of the Bareiss 
algorithm (6.1) can be written as follows: 



t 



(i+l) Ai+l) 



+2 



i+3 



1 



1 



L i+2 L i+3 
A-i) A-i) 
b i+2 L i+3 



A-i) 



(6.3) 
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where (4+2 > 4+2 > • • • ,4i ) are the last (n — i — 1) components of the first row of T(- f ), and 

(4+2)4+3' ■ ■ ■ ) 4^) are the last (n — i — 1) components of the first row of T«. 

Note that (6.3) has the same form as (3.13a)-(3.13b), and hence a connection between the 
Bareiss algorithm and algorithm FACTOR(T) is evident. That such a connection exists was 
observed by Sweet [22], and later by Delosme and Ipsen [11]. Sweet [22] related a step of the 
Bareiss algorithm (6.3) to a step of Bennett's downdating procedure [3]. Next, he derived the 
LU factorization of a Toeplitz matrix as a sequence of Bennett's downdating steps. Finally, he 
estimated the forward error in the decomposition using Fletcher and Powell's methodology [12]. 
This paper generalizes and presents new derivations of the results obtained in [22]. 



7 Numerical examples 

We adopt from [17] the following definitions of forward and backward stability. 

Definition 7.1: An algorithm for solving the equation (1.1) is forward stable if the computed 
solution x satisfies 

\\x — x\\ < ci(n)econd(T)||x|| , 

where cond(T) = ||T|| ||T _1 || is the condition number of T, and c\(n) may grow at most as fast 
as a polynomial in n, the dimension of the system. 

Definition 7.2: An algorithm for solving the equation (1.1) is backward stable if the com- 
puted solution x satisfies 

||Tx-6|| < c 2 (n)e||T||||x|| , 

where 02(11) may grow at most as fast as a polynomial in n, the dimension of the system. 

It is known that an algorithm (for solving a system of linear equations) is backward stable 
iff there exists a matrix AT such that 

(T + AT)x = b , 1 1 AT| I < c 3 (n)e||T|| , 

where c^(n) may grow at most as fast as a polynomial in n. 

Note that our definitions do not require the perturbation AT to be Toeplitz, even if the 
matrix T is Toeplitz. The case that AT is Toeplitz is discussed in [13, 24]. The reader is 
referred to [9, 14, 19] for a detailed treatment of roundoff analysis for general matrix algorithms. 

It is easy to see that backward stability implies forward stability, but not vice versa. This is 
manifested by the size of the residual vector. 

Cybenko [10] showed that the L\ norm of the inverse of a n x n symmetric positive definite 
Toeplitz matrix T n is bounded by 

-In , TT 1+ SUl! 

max-! ' ' 1 "' 11 ! 



'Wees* ft ' rea+smft)} - l|r " 1 ' 11 - n 1 



8111 t>i 

where {— shift}™", 1 are quantities called reflection coefficients. It is not difficult to pick the 
reflection coefficients in such a way that the corresponding Toeplitz matrix T n satisfies 

cond(T n ) wl/e. 

One possible way of constructing a Toeplitz matrix with given reflection coefficients {— sin ft}™^ 1 
is by tracing the elementary downdating steps backwards. 

An example of a symmetric positive definite Toeplitz matrix that can be made poorly con- 
ditioned by suitable choice of parameters is the Prolate matrix [21, 23], defined by 

f 2co if fc = 0, 

{ — otherwise, 
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where < oj < ^. For small to the eigenvalues of the Prolate matrix cluster around and 1. 

We performed numerical tests in which we solved systems of Toeplitz linear equations using 
variants of the Bareiss and Levinson algorithms, and (for comparison) the standard Cholesky 
method. The relative machine precision els 6 — _ 

10 16 . We varied the dimension of the 
system from 10 to 100, the condition number of the matrix from 1 to e _1 , the signs of reflection 
coefficients, and the right hand side so the magnitude of the norm of the solution vector varied 
from 1 to e _1 . In each test we monitored the errors in the decomposition, in the solution vector, 
and the size of the residual vector. 

Let xb and xl denote the solutions computed by the Bareiss and Levinson algorithms. Also, 
let t\b = Txb — b and rj_, = Txi — b. Then for the Bareiss algorithms we always observed that 
the scaled residual 

Iksll 
SB = 1 — iiFFir 

e||xB||||T|| 

was of order unity, as small as would be expected for a backward stable method. However, we 
were not able to find an example which would demonstrate the superiority of the Bareiss algo- 
rithm based on mixed downdating over the Bareiss algorithm based on hyperbolic downdating. 
In fact, the Bareiss algorithm based on hyperbolic downdating often gave slightly smaller errors 
than the Bareiss algorithm based on mixed downdating. In our experiments with Bareiss algo- 
rithms, neither the norm of the error matrix in the decomposition of T nor the residual error in 
the solution seemed to depend in any clear way on n, although a quadratic or cubic dependence 
would be expected from the worst-case error bounds of Theorems 5.1-5.2 and Corollary 5.1. 

For well conditioned systems the Bareiss and Levinson algorithms behaved similarly, and 
gave results comparable to results produced by a general stable method (the Cholesky method). 
Differences between the Bareiss and Levinson algorithms were noticeable only for very ill- 
conditioned systems and special right-hand side vectors. 

For the Levinson algorithm, when the matrix was very ill-conditioned and the norm of the 
solution vector was of order unity (that is, when the norm of the solution vector did not reflect 
the ill-conditioning of the matrix), we often observed that the scaled residual 

\\ r L\\ 

SL = 



e\\xL\\\\T\\ 

was as large as 10 5 . Varah [23] was the first to observe this behavior of the Levinson algorithm on 
the Prolate matrix. Higham and Pickering [16] used a search method proposed in [15] to generate 
Toeplitz matrices for which the Levinson algorithm yields large residual errors. However, the 
search never produced sl larger than 5 • 10 5 . It plausible that sl is a slowly increasing function 
of n and 1/e. 

Tables 7.1-7.3 show typical behavior of the Cholesky, Bareiss and Levinson algorithms for 
ill-conditioned Toeplitz systems of linear equations when the norm of the solution vectors is of 
order unity. The decomposition error was measured for the Cholesky and Bareiss algorithms by 
the quotient \\T — L • L T \\/(e • ||T||), where L was the computed factor of T. The solution error 
was measured by the quotient | \x comp — x\ |/||a;||, where x comp was the computed solution vector. 
Finally, the residual error was measured by the quotient \\T ■ x com p — feH/dl? 1 ) 





decomp. error 


soln. error 


resid. error 




Cholesky 
Bareiss (hyp) 
Bareiss (mixed) 
Levinson 


5.09 • KT 1 
3.45 • 10° 
2.73 • 10° 


7.67 • 1(T 3 

1.40 • 10" 2 

1.41 • 10° 
5.30 • 10° 


1.25 • 10° 
8.72 • 10" 1 
1.09 • 10° 
4.57- 10 3 
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Table 7.1: Prolate matrix, n = 21, w = 0.25, cond = 3.19 • 10 





decomp. error 


soln. error 


resid. error 




Cholesky 
Bareiss (hyp) 
Bareiss (mixed) 
Levinson 


1.72 • 10" 1 
2.91 • 10° 
3.63 • 10° 


6.84 • 10~ 2 
2.19 • 10" 1 
2.48 • 10" 1 
5.27- IO" 1 


3.11 • 10" 1 
1.15 • 10" 1 
2.47- 10" 1 
1.47- 10 5 





Table 7.2: Reflection coefficients | sin 9i\ of the same magnitude \K\ but 
alternating signs, \K\ = 0.8956680108101296, n = 41, cond = 8.5 • 10 15 





decomp. error 


soln. error 


resid. error 




Cholesky 
Bareiss (hyp) 
Bareiss (mixed) 
Levinson 


8.51 • 10" 1 
8.06 • 10° 
6.71 • 10° 


3.21 • 10~ 2 
1.13 • 10" 1 
1.16 • 10" 1 
2.60 • 10" 1 


4.28 • 10" 1 
2.28 • 10" 1 
3.20 • 10" 1 
1.06 • 10 5 





Table 7.3: Reflection coefficients | sin#j| of the same magnitude but 
alternating signs, \K\ = 0.9795872473975045, n = 92, cond = 2.77 • 10 



8 Conclusions 

This paper generalizes and presents new derivations of results obtained earlier by Sweet [22]. 
The bound in Corollary 5.1 for the case of mixed downdating is stronger than that given in [22]. 
The applicability of the Bareiss algorithms based on elementary downdating steps is extended to 
a class of matrices, satisfying Definition 4.2, which includes symmetric positive definite Toeplitz 
matrices. The approach via elementary downdating greatly simplifies roundoff error analysis. 
Lemmas 5.1 and 5.2 appear to be new. The stability of the Bareiss algorithms follows directly 
from these Lemmas and the results on the roundoff error analysis for elementary downdating 
steps given in [6]. 

The approach via downdating can be extended to the symmetric factorization of positive 
definite matrices of displacement rank k > 2 (satisfying additional conditions similar to those 
listed in Definition 4.2); see [18]. For matrices of displacement rank k the factorization algo- 
rithm uses elementary rank-A: downdating via hyperbolic Householder or mixed Householder 
reflections [8, 20]. 

We conclude by noting that the Bariess algorithms guarantee small residual errors in the 
sense of Definition 7.2, but the Levinson algorithm can yield residuals at least five orders of 
magnitude larger than those expected for a backward stable method. This result suggests that, 
if the Levinson algorithm is used in applications where the reflection coefficients are not known 
in advance to be positive, the residuals should be computed to see if they are acceptably small. 
This can be done in 0(n log n) arithmetic operations (using the FFT). 

It is an interesting open question whether the Levinson algorithm can give scaled residual 
errors which are arbitrarily large (for matrices which are numerically nonsingular) . A related 
question is whether the Levinson algorithm, for positive definite Toeplitz matrices T without a 
restriction on the reflection coefficients, is stable in the sense of Definitions 7.1 or 7.2. 
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